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A Symplectic Method to Generate Multivariate Normal Distributions 

C. Baumgarten 1 

Paul Scherrer Institute, Switzerland^ 
(Dated: 4 March 2013) 

The AMAS group at the Paul Scherrer Institute developed an object oriented library for high performance 
simulation of high intensity ion beam transport with space charge^. Such particle-in-cell (PIC) simulations 
require a method to generate multivariate particle distributions as starting conditions. 

In a preceeding publications it has been shown that the generators of symplectic transformations in two 
dimensions are a subset of the real Dirac matrices (RDMs) and that few symplectic transformations are 
required to transform a quadratic Hamiltonian into diagonal forir£ii. 

Here we argue that the use of RDMs is well suited for the generation of multivariate normal distributions 
with arbitrary covariances. A direct and simple argument supporting this claim is that this is the "natural" 
way how such distributions are formed. The transport of charged particle beams may serve as an example: 
An uncorrelated gaussian distribution of particles starting at some initial position of the accelerator is subject 
to linear deformations when passing through various beamline elements. These deformations can be described 
by symplectic transformations. 

Hence, if it is possible to derive the symplectic transformations that bring up these covariances, it is 
also possible to produce arbitrary multivariate normal distributions without Cholesky decomposition. The 
method allows the use of arbitrary uncoupled distributions. The functional form of the coupled multivariate 
distributions however depends in the general case on the type of the used random number generator. Only 
gaussian generators always yield gaussian multivariate distributions. 

PACS numbers: 45.20.Jj, 05.45.Xt, 41.85.-p, 02.50.-r 
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I. INTRODUCTION 

In Ref4 the author presented a so-called "decoupling" 
method that is based on the systematic use the real Dirac 
matrices (RDMs) in coupled linear optics. The RDMs 
are constructed from four pairwise anti-commuting basic 
matrices with the "metric tensor" = Diag(— 1, 1, 1, 1), 
formally written as: 



7m 7* +7^ 7m = 2 3Mf • 



(1) 



The remaining 12 RDMs are constructed as products of 
the basic matrices as described in the appendix. 

The use of the RDMs enables to derive a straightfor- 
ward method to transform transport matrices, force ma- 
trices ("symplices") and cr-matrices in such a way that 
the transformed variables are independent, i.e. decou- 
pled. 

The reverse is required to generate multivariate nor- 
mal distributions: A transformation that transforms lin- 
ear independent distributions of variables in such a way 
that a given covariance matrix is generated. The idea 
therefore is the following: Generate a set of independent 
normally distributed variables with given variances and 
apply the inverse of the decoupling transformation de- 
rived from the desired covariance matrix. This will cou- 
ple the "independent" variables in exactly the desired 
way. The presented scheme assumes an even number of 



variables since it is based on canonical pairs, i.e. position 
<7j and momentum pi - but it is always possible to ignore 
one of those variables. 

Since the method is based on pairs of canonical vari- 
ables, the decoupling scheme always treats two pairs of 
variables at a time, resulting in the use of 4 x 4-matrices. 
If more than four random variables are required, the de- 
coupling can be used iteratively in analogy to the Jacobi 
diagonalization scheme for symmetric matrices^. 



II. COUPLED LINEAR OPTICS 

In this section we give a brief summary of the major 
concept. Given the following Hamiltonian function 



(2) 



where A is a symmetric matrix and ip is a state-vector 
or "spinor" of the form ip = (qi,pi, q2,P2) T - The state 
vector hence contains two pairs of canonical variables. 
The equations of motion (EQOM) then have the familiar 
form 



A. = §H 
V = -M 



(3) 



or in vector notation: 
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where the force matrix F is given as F = 70 A. The ma- 
trix 70 is the symplectic unit matrix (sometimes labeled 
J or S) and is identified with the real Dirac matrix 70 
(see appendix) . We define the symmetric matrix of sec- 
ond moments tr containing the variances as diagonal and 
the covariances as off-diagonal elements. The matrix S 
is simply defined as the product of a with 70: 



S = a 70 . 



(5) 



Both matrices, F and S, fulfill the following equation 
(using 7^ = -70 and 7^ = -1): 



70 F 70 . 



(6) 



Matrices that obey Eq. [6] have been named symplices, 
but they are also called "infinitesimally symplectic" or 
"Hamiltonian" matrices 5 . Symplices allow superposition, 
i.e. any sum of symplices is a symplex, but only the 
product of anti-commuting symplices is a symple>s^. 

Any real- valued 4 x 4-matrix M can be written as a 
linear combination of real Dirac matrices (RDM): 



M 



15 

^ m k 7 fe 

k=0 



(7) 



The RDM-coefficients nik can be computed from the ma- 
trix M by: 



Tr( 7fc 2 ) 
32 



Tr(M 7fc + 7fe M), 



(8) 



where Tr(X) is the trace of X. 

Hence the RDMs form a complete system of all real 
4 x 4-matrices, but only ten RDMs fulfill Eq. [6] and are 
therefore symplices: The basic matrices 70,... ,73 and 
the six "bi- vectors" , i.e. the six possible products of two 
basic matrices. The symplices are the generators of sym- 
plectic transformations, i.e. the generators of the sym- 
plectic group. 

As well-known, the Jacobi matrix of a canonical trans- 
formation is symplectic, i.e. it fulfills the following equa- 
tion^: 



M 70 M J = 70 . 
The EQOM have the general solution 

^(t) = M(t,i<,)V(to), 



(9) 



(10) 



where M is a symplectic transfer matrix that is in case 
of constant forces given by 



M(Mo) =exp(F (t-t )). 



(11) 



Given now an (initial) set of N normally distributed 
uncorrelated random variables ipi, then the tr-matrix of 
these variables is given by 



N-l 
i=0 



(12) 



where the superscript "T" indicates the transpose, then 
the distribution at time t is given by: 



N-l 



at = Jq 2 M^V-f M T = M(7oM T . (13) 



Hence with Eqn. §5§ and © one has: 



St = -M f j 7o 2 M T 7o 
= MSqM- 1 . 



(14) 



That is - the transformation of S is a similarity- 
transformation with a symplectic transformation matrix. 
The reverse transformation obviously is 



S = M-'SfM. 



(15) 



Now we refer to the structural identity of the matrix S 
with the force matrix F. Both are symplices and since 
a transformation that decouples F has been shown to 
diagonalize the matrix A of the Hamiltonian^ii, it is clear 
that the same method can be used to diagonalize a. The 
reverse of this transformation then generates the desired 
distribution from an initially uncorrelated a. 

Instead of a Cholesky-decomposition we may therefore 
use a symplectic similarity-transformation to generate 
the correlated distribution from an initially uncorrelated 
distribution. In the context of charged particle optics, 
the algorithm delivers even more useful information: the 
transformation matrix M _1 is the transport matrix that 
is required to generate an uncorrelated beam. 



III. SYMPLECTIC TRANSFORMATIONS AND THE 
ALGORITHM 

The general form of a symplectic transformation ma- 
trix R/, is that of a matrix exponential of a symplex 75 
multiplied by a parameter e representing cither the angle 
or the "rapidity" : 



R*(e) 
-i, 

b 



exp( 7b §) = lc + j b s 
cxp(- 7fc |) =lc-7 6 s, 



(16) 



where 



cos(e/2) for 

cosh(e/2) for 

sin(e/2) for 

sinh(e/2) for 



if, 



1 
1 
1 
1 



(17) 



Transformations with 7 ^ = — 1 are orthogonal transfor- 
mations, i.e. rotations, while those with 7^ = 1 are 
boosts. 

The matrix S then is transformed according to: 



RSR 



(18) 
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The decoupling requires a sequence of transformations, 
so that the RDM-coefficients of S have to be recomputed 
after each step. 

Eqn. [5] may be used to compute the RDM-coefficients 
Sk of the matrix S 



9 
i=0 



(19) 



Numerically it is faster to analyze directly the composi- 
tion. For the choice of RDMs used in RefA 4 the RDM- 
coefficients of S as a function of a are given by: 



•so 


= (on + 


0"22 + CT33 + 


(744 )/4 


Si 


= (-o-ii 


+ CT 2 2 + O33 


- 0-44 )/4 


S2 


= (o"13 - 


ct 24 )/2 




S3 


= (0"12 + 


t 34 )/2 




S4 


= (0"12 - 


T34)/2 




S5 


= -(014 


+ ct 23 )/2 


CT 4 4)/4 


S6 


= (o"n - 


0"22 + C33 — 


S7 


= (0"13 + 


<724)/2 


CT 44 )/4 


S8 


= (on + 


0"22 — C33 — 


S9 


= (C14 - 


a 23 )/2 





(20) 



Now we use the following abbreviation using the notation 
of 3-dimensional vector algebra: 



S 
P 
E 
B 



and furthermore: 



M r 


= EB 


f — 


M g 


= BP 


9 = 


M b 


= EP 


b = 



so 

(si,s 2 ,s 3 ) T 
(s 4 ,s 5 ,s 6 ) T 

(S7, Sg, S 9 ) T . 



SP- 
SS- 



(21) 



B x E 
P x B 
ExP 



(22) 



The decoupling is done by a sequence of maximal six 
symplectic transformations 4 . A transformation with e = 
can be omitted. After each transformation, the RDM- 
coefficients Sk have to be updated and Eqns. (|2~T1) and 
(1221) have to be re-evaluated: 

1. Ro(e) with e — arctan(^-). 



2. Kr(e) with e = arctan(^) 

3. Rg(£) with e 

4. R 2 (e) with e 

5. Ro(e) with e = h arctan( 

6. Rs(e) with e 



— arctan (jp). 
artanh(^t). 



2M b 



E 2 -P 2 ' 



— arctan (jr). 



Given an initial covariance matrix ao; the sequence of 
computation therefore is: 



1. Compute the RDM-coefficients Sk according to 
Eqn. (|2U)) and the quantities defined in Eqns. (|2"Tj) 
and (1221). 



2. Compute the first (or next, resp.) transformation 
matrix R. 

3. Compute the product of the transformation matri- 
ces (and of the inverse) M„+i = R n +i R«- 

4. Apply the first (or next, resp.) transformation 
S?i+i = R S„ R 

5. Compute a n+ i = -S n+ i7 . 

6. Continue with next transformation at step 1). 

The six iterations yield the desired diagonal matrix <76 
and the matrices M 6 and its inverse, so that 



M 6 S M 6 



(23) 



(24) 



The diagonal elements of eg are the variances of the un- 
coupled gaussian distribution. Given ipi is the i-th un- 
coupled random state vector, then M^ 1 tpi is the corre- 
sponding state vector with the multivariate normal dis- 
tribution. 



IV. EXAMPLE 

Consider for instance the (arbitrary) matrix of second 
moments ctq 



5.8269 
-0.0303 

0.2292 

0.0000 
-0.0960 

1.4897 



-0.0303 
0.8851 
0.0000 

-0.0311 
1.8053 

-0.0015 



0.2292 
0.0000 
3.6058 
-0.0235 
0.0000 
0.0000 



0.0000 
-0.0311 
-0.0235 
0.6844 
0.0000 
0.0000 



-0.0960 
1.8053 
0.0000 
0.0000 
7.0607 

-0.0224 



The diagonal matrix ae is computed to be 



1.4897 \ 
-0.0015 

0.0000 

0.0000 
-0.0224 

0.7304 , 
(25) 



/ 1.9982 




V 





1.2984 










1.1116 










2.2124 










1.4029 










1.6637 / 



(26) 



Now 10 5 random vectors have been generated with a 
Gaussian random number generator of unit variance. 
The vector elements have been scaled with corresponding 
variances, given by the root of the diagonal elements of 
cr 5 and then been multiplied (or transformed) with M _1 
given by 



-0.1727 

-0.0371 

0.9474 

-0.0573 

-0.1641 

-0.3455 



-0.0330 
0.3392 
0.1485 
0.5025 
1.6387 

-0.0555 



0.0081 
0.7051 
0.0008 
-0.0072 
0.3732 
0.0042 



-1.6049 
0.0093 

-0.0613 
0.0001 
0.0202 

-0.3515 



-0.0725 
0.3402 

-0.3429 

-0.4733 
1.4755 

-0.1204 



0.1893 \ 
0.1034 
0.9838 
-0.1464 
0.4320 
0.3416 , 
(27) 
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FIG. 1. Top: Correlations between several variables of 
the multivariate normal distribution. Bottom: The result- 
ing probability distributions for the individual variables are 
again gaussian. 



Then the covariance matrix of the produced random 
vectors was evaluated. The result is: 



/ 5.7946 -0.0247 0.2343 -0.0060 -0.1036 1.4825 \ 
-0.0247 0.8917 0.0023 -0.0306 1.8200 0.0034 
0.2343 0.0023 3.5910 -0.0299 -0.0158 0.0033 
-0.0060 -0.0306 -0.0299 0.6849 -0.0000 -0.0012 
-0.1036 1.8200 -0.0158 -0.0000 7.0928 -0.0148 
\ 1.4825 0.0034 0.0033 -0.0012 -0.0148 0.7294 / 

(28) 

Fig. [T] shows some of the distributions as examples. 

The same procedure can be done with any initial prob- 
ability distribution and the algorithm will produce the 
desired second moments. But the functional form of the 
resulting distributions of the transformed variables will 
only be similar to the initial distribution in the Gaussian 
case. Fig. [5] shows the results for the same covariance 
matrix if the decoupled variables have a uniform proba- 
bility distribution, but same variances. The covariance 
matrix is correctly reproduced. 



V. CONCLUSION 

The method of symplectic decoupling of linearily cou- 
pled variables has been applied to the problem of mul- 
tivariate random distributions. It has been shown that 
the use of sympleptic algebra has severe advantages: The 
same methods can be applied to solve a variety of prob- 
lems. The presented algorithm is especially interesting 
for the generation of starting conditions of particle track- 
ing codes like - for example - OPALi^. 

In cases where the decoupled process is known to have 
a non-Gaussian probability distribution and if the trans- 
port matrix M of a linear transport system is known, it 
should be possible to derive unknown parameters of the 
initial distribution by comparison with the computed ex- 
pected distribution. Fig. [5] shows that a fiat distribution 
yields a clear "signature" . 



FIG. 2. Top: Correlations between several variables of the 
multivariate flat distribution. Bottom: The resulting prob- 
ability distributions for the individual variables strongly de- 
pend on the correlations. The stronger the correlations, the 
more gaussian the distribution will be (central limiting theo- 
rem). 
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Appendix A: The 7-Matrices 

The real Dirac matrices used throughout this paper 
are: 



70 



72 

714 
74 
75 

76 

710 
712 




70717273; 
7071; 
7072; 
7073; 

714 70 = 7l 72 73 
714 72 = 70 73 71 



71 



73 

715 

77 

78 

79 
711 
713 





-1 




-1 






. 



1 




-1 





714 70 71 = 72 73 
714 70 72 = 73 71 
714 7o 73 = 71 72 
714 71 = 7o 72 73 
714 73 = 7o 7i 72 

(Al) 
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